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TWO-EQUATION,  DEPTH- INTEGRATED  TURBULENCE  CLOSURE 
FOR  MODELING  GEOMETRY-DOMINATED  FLOWS 


Introduction 


1.  When  applying  depth-integrated  hydrodynamic  and  dispersion 
models  to  investigate  large-scale  water  quality  problems  in  coastal 
environments,  it  is  usually  adequate  to  adopt  the  simplest  of  eddy 
viscosity/diffusivity  closure  hypotheses  in  which  the  depth-integrated 
eddy  viscosity/diffusivity  coefficient  is  either  assumed  to  be  constant 
or  related  to  the  water  depth  and  shear  velocity.  However,  when  local 
geometry-dependent  flow  phenomena  such  as  separation  and  recirculation 
are  deemed  important  to  the  water  quality  investigation,  it  is  then 
necessary  to  adopt  a  more  sophisticated  level  of  closure  for  the  depth- 
integrated  equations  of  motion  and  constituent  transport  equations. 

2.  It  is  the  intent  of  this  discussion  to  (a)  describe  the  clo¬ 
sure  problem  associated  with  the  use  of  the  depth-integrated  equations 
of  motion,  (b)  provide  a  detailed  description  of  the  two-equation  (k-e) 
turbulence  closure  model,  and  (c)  present  the  results  of  steady-state 
test  simulations  utilizing  k-e  turbulence  closure. 

3.  In  addition,  the  limitations  of  the  model  developed  by  the 
author  are  discussed  with  specific  recommendations  for  future 
improvement . 


Turbulence  Closure 


Depth-integrated  equations 
of  motion  and  closure  problem 

4.  Under  the  assumption  of  a  homogeneous,  incompressible,  viscous 
flow  characterized  by  a  hydrostatic  pressure  distribution,  with  wind  and 
Coriolis  forces  neglected,  the  depth-integrated  equations  of  motion  are 
written: 


Conservation  of  Mass 
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where 

m,n  =1,2,  and  repeated  indices  require  summation 
t  =  time 

V  =  two-dimensional  depth-averaged  velocity  vector  (U,  V) 
h  =  water  depth 

X  =  coordinate  directions  (x,  y) 
m 

g  =  acceleration  due  to  gravity 

*  components  of  the  bottom  shear  stress  per  unit  mass 
«  channel  bottom  elevation  above  an  arbitrary  datum 

T  =  components  of  the  depth-integrated  effective  stress  tensor 
per  unit  mass 

5.  The  depth-integrated  effective  shear  stress  tensor  as  defined 
by  Kuipers  and  Vreugdenhil  (1973)  and  Flokstra  (1977)  contains  the 
(I)  viscous  stresses,  (II)  the  turbulent  Reynold's  stresses,  and 
(III)  the  momentum  dispersion  terms  which  arise  from  depth  integrating 
the  nonlinear  convective  acceleration  terms  in  the  equations  of  motion. 
Specifically,  the  depth-integrated  effective  stress  tensor  per  unit  mass 
is  written: 
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where 
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v  =  kinematic  viscosity 

=  time-averaged  velocity  components  (u,  v) 

v^  =  horizontal  turbulent  velocity  fluctuations 

6.  The  closure  problem  associated  with  the  use  of  the  depth- 
integrated  equations  of  motion  results  from  the  need  to  parameterize 
terms  II  and  III  in  the  effective  stress  tensor. 

7.  The  contribution  of  the  viscous  stresses  (term  I)  can  be 
neglected  simply  because  its  effects  are  transparent  to  the  computation 
at  the  scales  of  motion  modeled  in  a  fully  turbulent  hydraulic  problem. 

8.  In  the  recent  literature,  the  treatment  of  Reynolds  stress 
closure  (term  II)  has  ranged  from  neglecting  the  terms  to  the  concept  of 
"large  eddy"  simulations  (Leonard  1974)  where,  by  means  of  spatial 
filtering  of  the  equations  of  motions,  only  the  small-scale  or  "subgrid- 
scale"  Reynolds  stresses  need  be  modeled.  Within  these  two  extremes  lie 
a  number  of  alternative  closure  schemes  which  exhibit  a  wide  variation 
in  complexity  (Reynolds  1976). 

9.  One  closure  technique  that  has  enjoyed  considerable  success  in 
the  simulation  of  a  variety  of  turbulent  flows  is  the  two  equation  (k-e) 
turbulence  model  described  by  Launder  and  Spalding  (1974)  and  Rodi 
(1980).  Applications  of  a  depth-integrated  version  of  the  k-t  turbu¬ 
lence  model  have  been  presented  by  Rastogi  and  Rodi  (1978)  and  McQuirk 
and  Rodi  (1978) .  A  detailed  discussion  of  the  three-dimensional  and 
depth-integrated  version  of  the  k-e  closure  model  is  the  subject  of 
the  next  section. 

10.  Parameterization  of  the  momentum  dispersion  terms  (III)  is 
straightforward  if  one  has  an  a  priori  knowledge  of  the  vertical  dis¬ 
tribution  of  the  horizontal  velocity  components.  Unfortunately,  theo¬ 
retical  velocity  distributions  are  available  for  only  the  simplest  of 
flows  such  as  flow  in  a  wide  channel  or  long  circular  channel  bend. 
Discussions  of  the  importance  of  the  momentum  dispersion  terms  and 
approximate  closure  schemes  are  presented  by  Flokstra  (1977),  Abbott  and 
Rasmussen  (1977),  and  Lean  and  Weare  (1979).  A  common  feature  of  all  of 
the  existing  closure  schemes  for  momentum  dispersion  is  that  the  magni¬ 
tude  of  the  components  of  momentum  dispersion  that  can  be  important  is 


directly  proportional  to  the  ratio  of  the  depth  of  flow  to  the  radius  of 
curvature  of  the  depth-mean  streamlines.  Consequently,  if  one  restricts 
their  attention  to  flows  that  are  much  wider  than  deep,  which  is  usually 
the  case  in  estuarine  and  other  coastal  flows,  then  the  need  to  address 
momentum  dispersion  is  obviated. 

Depth-averaged  turbulent 
Reynold's  stress  closure 

11.  The  k-e  turbulence  closure  model  presented  by  Launder  and 
Spalding  (1974)  is  based  on  the  Boussinesq  eddy  viscosity  hypothesis 
(Hinze  1959),  which  assumes  that  the  turbulent  Reynold's  stresses  are 
proportional  to  the  mean  strain  rates.  Using  three-dimensional  tensor 
notation,  the  turbulent  Reynold's  stresses  are  written: 


where 

i>j  =  1>  2,  3,  and  repeated  indices  require  summation 
k  =  ~2  vivi’  t*le  turbulent  kinetic  energy  per  unit  mass 
6^  =  Kronecker  delta 

Unlike  the  molecular  viscosity  v  ,  the  turbulent  eddy  viscosity  vt 
is  flow-dependent  and  can  vary  both  in  space  and  time.  An  approxima¬ 
tion  for  the  distribution  of  the  turbulent  eddy  viscosity  is  obtained  by 
assuming  that  it  is  proportional  to  the  product  of  the  characteristic 
velocity  and  length  scale  of  turbulence,  namely: 

,  1/2 . 

vt  "  k  i  (5) 

where  i  equals  the  macroscale  of  turbulence  (a  measure  of  the  size  of 
the  energy  containing  eddies) .  An  inviscid  estimate  of  the  energy 
dissipation  rate  per  unit  mass  e  is  obtained  when  one  assumes  that  the 
amount  of  energy  dissipated  at  the  small  scales  of  turbulence  equals  the 
rate  of  supply  at  the  large  scales. 
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12.  Again  utilizing  the  characteristic  velocity  and  length  scales 
of  turbulence,  dimensional  considerations  require  that  (Tennekes  and 
Luraley  1972) 
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k3/2 
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(6) 


Substitution  of  Equation  6  into  Equation  5  yields  a  functional  relation¬ 
ship  for  the  turbulent  eddy  viscosity  in  terms  of  the  kinetic  energy  of 
turbulence  k  ,  and  its  rate  of  dissipation  e  ,  specifically: 


v 

t 


v  e 


(7) 


where  is  an  empirical  coefficient. 

13.  In  principle,  solution  of  the  exact  transport  equations  for 
the  turbulence  kinetic  energy  k  ,  and  its  rate  of  dissipation  e  , 
enables  one  to  completely  specify  the  temporal  and  spatial  distribution 
of  the  turbulent  eddy  viscosity.  However,  construction  of  the  transport 
equations  for  k  and  e  results  in  an  additional  closure  problem. 
Consider  first  the  turbulence  kinetic  energy  equation  (Hinze  1959)  : 


(III)  (IV) 


where  P*  denotes  turbulent  pressure  fluctuations,  P  is  the  fluid 
density,  and  the  over-bar  represents  a  time  average.  The  closure  prob¬ 
lem  results  from  the  presence  of  the  unknown  pressure  and  velocity 
fluctuation  correlations  in  term  I.  Physically,  term  I  represents  the 
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convective  diffusion  of  the  total  turbulence  mechanical  energy  per  unit 
mass  by  turbulence.  This  term  acts  as  a  redistribution  mechanism,  which 
suggests  the  use  of  a  gradient  diffusion  model  (Rodi  1980),  or: 


i  \ P  /  ak  3xi 


(9) 


where  is  an  empirical  constant.  Equation  4  may  be  substituted 

directly  for  the  turbulent  Reynold's  stresses  in  the  turbulence  produc¬ 
tion  term  (II).  Term  IV  is  by  definition  the  energy  dissipation  rate 
per  unit  mass,  e  .  Finally,  term  III  represents  the  work  done  by  the 
viscous  shear  stresses.  For  high  Reynold's  number  flows,  this  term  is 
small  and  can  be  neglected  (Hanjalic  and  Launder  1972).  Making  the 
appropriate  substitutions,  the  three-dimensional  model  equation  for  the 
turbulence  kinetic  energy  is  written: 


3(v.k) 


3x 
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14.  The  exact  transport  equation  for  the  energy  dissipation  rate 
per  unit  mass,  for  large  Reynold's  numbers,  reads  (Harlow  and  Nakayama 
1968)  : 


(IV)  (V) 
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t 

I 

where 

r,s  =  1,  2,  3,  and  repeated  indices  require  summation 

s'  =  turbulent  fluctuations  of  the  energy  dissipation  rate  per 
unit  mass 

15.  The  closure  approximations  for  Equation  11  were  first  pre¬ 
sented  by  Hanjalic  and  Launder  (1972).  Their  approach  was  to  param¬ 
eterize  I-III  in  terms  of  the  Reynold's  stresses,  mean  strain  rate, 
turbulence  kinetic  energy,  and  its  rate  of  dissipation  per  unit  mass, 
and  to  neglect  term  V  on  the  basis  of  being  small.  Term  I  represents 
the  production  mechanisms  for  e  and  was  approximated  accordingly 

Term  I  =  -C,  (£)  (v*v;>  —  (12) 

\  I  r 

Terms  II  and  III  were  grouped  together  and  parameterized  as  follows: 

2 

II  +  III  -  C2  ~  (13) 

16.  The  argument  used  to  support  Equation  13  was  that  the  sum  of 
term  II,  which  represents  the  generation  rate  of  vorticity  fluctuations 
due  to  the  self-stretching  mechanism,  and  term  III,  which  represents  the 
viscous  decay  of  dissipation,  should  be  controlled  by  the  dynamics  of 
the  energy  cascade.  Subsequently,  if  the  Reynold's  number  is  suffi¬ 
ciently  large  to  allow  the  existence  of  an  inertial  subrange,  an  invis- 
cid  estimate  based  on  dimensional  considerations  should  be  appropriate. 
Term  IV  represents  the  turbulent  diffusion  of  e  ,  which  clearly 
suggests  a  closure  of  the  form: 

v 

Term  IV  =  —  —  (14) 

a  x . 

e  j 

where  is  an  empirical  constant.  Collecting  the  various  approxima¬ 

tions  yields  the  complete  model  equation  for  the  energy  dissipation  rate 
per  unit  mass: 


9 
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17.  Estimates  for  the  empirical  constants  found  in  Equations  10 
and  15  were  obtained  by  applying  the  model  equations  to  simple  turbulent 
flows  for  which  experimental  data  were  available.  For  example,  in  a 
local  equilibrium,  two-dimensional  boundary  layer,  the  production  of 
turbulence  energy  is  balanced  by  dissipation,  and  Equation  10  reduces 


— i — — 

U  V  -r—  =  C 
3y 


A  consistent  closure  approximation  for  the  Reynold's  stress  is: 


— ; — r  3u 

u'v  =  v  — 


which  yields 


— i — r2 

U  V  *  VtE 


Now,  by  definition, 


v  =  C  — 
t  VC 


Thus,  substitution  of  Equation  18  into  Equation  19  results  in 


,1/2  _  u'v' 


18.  Experimental' data  extracted  from  the  work  of  Laufer  (1951) 
suggest  that  u’v’/k  varies  from  0.22  to  0.3,  which  gives  a  range  of 
from  0.05  to  0.09.  The  value  of  C ^  was  found  to  lie  between  1.9 
and  2.0  (Hanjalic  and  Launder  1972;  Launder  and  Spalding  1974)  when 
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computed  using  measured  decay  rates  of  the  turbulent  kinetic  energy 
behind  a  grid  (Townsend  1956) .  The  constant  was  obtained  by  exam¬ 

ining  the  form  of  the  dissipation  equation  in  the  constant  shear  stress 
region  near  the  wall.  Here,  convection  is  negligible  with  production 
and  dissipation  approximately  in  balance,  thus 


e 


<y 


(21) 


where 

U*  =  shear  velocity 
k  =  von  Kamans  constant 
y  =  distance  from  the  wall 

Substituting  Equation  21  into  Equation  15  with  some  simplification 
yields 
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In  the  near  wall  region  (Townsend  1956) 
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thus 
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where  by  definition 
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(24) 
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and  therefore 


which  specifies  the  value  of  when  ,  C2  *  CTk  ’  aru*  °t 

are  known.  The  constants  and  a ^  ,  which  are  similar  to  turbulent 

Prandtl  or  Schmidt  Numbers,  were  assumed  to  be  of  order  unity  and  deter¬ 
mined  along  with  the  recommended  values  of  the  other  constants  via 
computer  optimization.  This  was  done  by  adjusting  the  values  of  the 
various  constants  until  reasonable  agreement  between  computed  and  experi¬ 
mental  results  was  obtained.  The  values  recommended  by  Launder  and 
Spalding  (1974)  are  as  follows: 


Cv  -  0.09 
Cj  -  1.44 
C2  -  1.92 

ak  "  1*° 
a  -  1.3 

E 

With  the  appropriate  specifications  of  boundary  and  initial  conditions. 
Equations  7,  10,  and  15  and  the  Launder  and  Spalding  values  presented 
above  constitute  the  complete  three-dimensional  k-E  turbulence  closure 
model. 

19.  To  be  of  use  in  approximating  the  depth-averaged  Reynold's 
stress  in  Equation  7,  it  is  necessary  to  cast  the  three-dimensional 
k-E  model  into  a  depth-integrated  form.  Realizing  that  turbulence  is 
inherently  three-dimensional,  depth  integration  of  the  transport  equa¬ 
tions  for  the  turbulence  kinetic  energy  and  its  rate  of  dissipation 
cannot  be  strictly  performed.  However,  Rastogi  and  Rodi  (1978)  suggest 
model  equations  for  the  depth-averaged  turbulence  energy  k  ,  and  its 
rate  of  dissipation  e  ,  can,  in  fact,  be  constructed  if  additional 
source  terms  are  added  to  account  for  mechanisms  originating  from  non¬ 
uniformity  of  the  flow  over  the  vertical  dimension.  Furthermore,  they 
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suggest  that  the  resulting  turbulent  viscosity  v  should  be  inter¬ 
preted  such  that  when  multiplied  by  the  depth-averaged  strain  rate  will 
yield  the  depth-averaged  turbulent  Reynold's  stress.  By  analogy,  the 

depth-averaged  Reynold's  stress  tensor  I  is  written 

inn 


.  3(Vh)  3(V  h)1  „ 

Z  “  ~  -3  +  — -,v~  ■  -  4  kh6 

mn  t  3X  3X  3  mn 

n  m 


where 


v  =  C  — 
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20.  The  resulting  model  equations  for  the  depth-averaged  value  of 

«  A 

the  turbulence  energy  k  and  its  rate  of  dissipation  e  are  written  as 
follows: 


,  3(v  kh) 

3(kh)  .  m 

3t  3X 

m 


=  kmoI  +  p 

3X  t  3X  h 

m  m  , 


+  P.  -  eh 
k 


3  <eh)  +  „  JL  lL  iiilll  +£(Cp 

3t  3X  3X  a  3X  .*  l  l  h 

m  m  e  m  k 


C2eh)  +  Pe 


where 


v.  3 (V  h)  3(V  h)  3(V  h) 
p  „  _L  .  m  +  S  m 

h  h  3X  3X  3X 

n  m  .  n 


21.  The  source  terms  P^  and  P£  account  for  the  production 
mechanism  resulting  from  the  presence  of  a  vertical  boundary  layer.  The 
form  of  these  production  terms  may  be  obtained  by  considering  the 


i 


central  portion  of  a  unidirectional,  uniform  flow  in  a  wide  open  chan¬ 
nel.  For  this  flow,  the  balance  equations  for  k  and  e  reduce  to 


P 


k 


th 


(32) 


and 


£2h 


(33) 


To  a  good  approximation,  the  total  turbulence  energy  production  over  the 
vertical  is  written  (Townsend  1956)  : 


Pk  =  U*U  (34) 

but,  by  definition, 

U  l  =  cl)2  (35) 

where  c  is  a  nondimensional  friction  coefficient;  therefore, 

Pk  =  cU3  (36) 


A  similar  relation  may  be  obtained  for  the  dissipation  source  P^  by 
recalling  that  from  Equation  28 


which  may  be  substituted  into  Equation  33  to  yield 


P 


„  „l/2 '3/ 2. 

u  C  e  h 
_£v _ 


(37) 


(38) 
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Noting  that 


one  defines  D  to  be  a  longitudinal  dispersion  coefficient  for  an 
infinitely  wide  open  channel,  its  value  is  approximately  5.9.  Between 
these  two  extremes,  an  entire  spectrum  of  values  may  be  obtained  if  D 
is  considered  to  be  a  transverse  mixing  coefficient.  Depending  on  the 
cross-sectional  shape  and  the  longitudinal  curvature  of  the  channel 
investigated,  a  compilation  of  numerous  experiments  yields  values  rang¬ 
ing  from  0.1  to  1.0  (Fischer  et  al.  1979).  Nonetheless^  computational 
experience  has  shown  (Chapman  1982)  that  a  value  of  unity  for  the  non- 
dimensional  dispersion  coefficient  yields  satisfactory  results  in  the 
simulation  of  a  depth-integrated  wall  boundary  layer. 

Results  of  Steady-State  Test  Simulations 

24.  The  model  problem  chosen  for  the  test  simulations  was  flow  in 
a  wide,  shallow,  rectangular  channel  with  an  abrupt,  symmetric  expansion 
in  width  (Figure  1) .  The  reasons  for  adopting  this  test  problem  are 


essentially  twofold.  First,  a  range  of  reliable  nondimensional  reattach¬ 
ment  lengths  xr/wi  are  available  in  the  form  of  review  papers  (Kim, 
Kline,  and  Johnston  1978;  Durst  and  Tropea  1981;  Eaton  and  Johnston, 
1981).  Although  the  value  of  the  measured  nondimensional  reattachment 
varies  considerably  from  one  experiment  to  the  next,  it  was  first 
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pointed  out  by  de  Brederode  and  Bradshaw  (1972)  that  much  of  the  varia¬ 
tion  in  reattachment  length  measurements  could  be  attributed  in  differ¬ 
ences  in  the  aspect  ratios  of  the  test  sections.  The  observation  of 
de  Brederode  and  Bradshaw  is  well  illustrated  in  Figure  2,  a  plot  of 


h/W„ 


Figure  2.  Experimental  measurements  of  nondimensional  reattachment 
lengths  versus  the  inlet  channel  aspect  ratio 

published  measurements  of  reattachment  lengths  for  symmetric  channel 
expansions  as  a  function  of  the  upstream  section  aspect  ratio,  h/WQ  . 

In  this  figure  the  trend  of  decreasing  reattachment  length  with  decreas¬ 
ing  aspect  ratio  is  clearly  seen. 

25.  The  second  reason  for  choosing  the  channel  expansion  test 
problem  is  that  previous  attempts  at  simulating  this  flow,  using  simple 
eddy  viscosity  models,  predicted  reattachment  lengths  that  were  three  to 
four  times  too  short  (Abbott  and  Rasmussen  1977;  Ponce  and  Yabusaki 
1981). 

26.  The  details  of  the  test  simulations  are  not  presented  herein; 
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however,  they  may  be  found  in  Chapman  (1982)  or  Chapman  and  Kuo  (1982). 
The  results  of  the  initial  simulation  applying  the  standard  depth- 
integrated  k-e  turbulence  closure  model,  as  described  in  the  previous 
section,  are  presented  in  Figure  3,  a  vector  plot  of  the  depth-integrated 


I _ 1  RlNLtT  NOTE:  GRID  POINTS  ARE  LOCATED  AT  THE 

0  1  MIDPOINT  OF  THE  VELOCITY  VECTORS 


Figure  3.  Depth-averaged  velocity  field  for  the  standard 
(k-e)  turbulence  closure  simulation 

velocity  field.  In  this  simulation,  the  inlet  section  aspect  ratio  is 
about  0.1,  which  suggests  that  the  nondimensional  reattachment  length 
Xr^l  be  about  A. 5  to  5.0.  However,  it  is  seen  in  Figure  2  that 

the  predicted  reattachment  length  is  only  about  3.2,  which  corresponds 
to  an  error  of  about  30  percent.  The  poor  agreement  between  model 
prediction  and  experimental  measurements  is  directly  attributable  to 
dependence  or  the  coefficient  C  to  the  degree  of  curvature  of  the 
depth-mean  streamlines  (Leschziner  and  Rodi  1981) .  In  an  attempt  to 
improve  upon  the  model  predictions,  an  ad  hoc  approximation  to  the 
streamline  curvature  modification  of  Leschziner  and  Rodi  (1981)  was 
employed.  Specifically,  the  value  of  the  coefficient  C  was  decreased 
to  0.03  at  all  grid  points  in  the  region  behind  the  step  for  a  distance 
of  6W^  downstream. 

27.  The  results  of  the  curvature-corrected  simulation  are  pre¬ 
sented  in  Figure  A  in  which  a  significant  increase  in  the  length  of  the 
predicted  reattachment  length  is  seen.  The  value  of  the  predicted 
reattachment  length  is  about  x^W^  -  A. 6  ,  which  agrees  well  with  the 
experimental  measurements  depicted  in  Figure  2. 
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NOTE:  GRID  POINTS  ARE  LOCATED  AT  THE 

MIDPOINT  OF  THE  VELOCITY  VECTORS 


Figure  4.  Depth-averaged  velocity  field  for  curvature-corrected 
(k-e)  turbulence  closure  simulation 


Model  Limitations  and  Recommendations  for  Future  Work 


28.  The  failure  of  the  standard  depth- integrated  k-e  turbulence 
closure  formulation  to  predict  correct  reattachment  lengths,  and  the 
marked  improvement  that  was  realized  using  an  ad  hoc  streamline  curva¬ 
ture  correction,  suggests  that  the  complete  streamline  correction  formu¬ 
lation  of  Leschziner  and  Rodi  should  be  incorporated  into  the  model. 

29.  Secondly,  the  present  model  uses  a  simple  explicit  temporal 
uptake  procedure  to  iterate  the  solution  to  steady-state.  In  order  to 
apply  the  model  to  slowly  varying  environmental  flows  (i.e.,  tidal 
flows),  it  will  be  necessary  to  implement  a  version  in  which  the  water 
surface  elevation  is  computed  implicitly,  which  removes  the  overly 
restructive  gravity  wave  propagation  speed  stability  criteria  (Johnson 
1980) . 

30.  Thirdly,  if  one  desires  to  address  rapidly  varying  transient 
phenomena,  then  the  existing  spatially  third-order  numerical  technique 
QUICK  (Quadratic  Upstream  Interpolation  for  Convective  Kinematics)  of 
Leonard  (1979)  must  be  modified  in  such  a  way  that  it  will  also  be 
temporally  third  order.  The  resulting  algorithm,  which  is  called 
QUICKEST  (Quadratic  Upstream  Interpolation  for  Convective  Kinematics 
with  Estimated  Streaming  Terms) ,  was  presented  in  one  dimension  by  its 


■s 
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originator  Leonard  (1979)  and  extended  to  a  general  two-dimensional  form 
by  Hall  and  Chapman  (1982). 

31.  Finally,  at  the  present  time,  the  locations  of  boundaries  are 
programmed  into  the  mainline  of  the  code,  which  requires  reprogramming 
with  every  new  geometry  configuration.  Consequently,  it  will  be  neces¬ 
sary  to  incorporate  a  flag  system  that  will  enable  any  boundary  config¬ 
uration  to  be  generated  via  input  data. 
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